Same as Dataset 2f but sample size of 50 and number of taxa to 100.
All different interaction matrices
1. Load Library
library(tidyverse)
library(igraph)
library(NBZIMM)
library(SpiecEasi)
library(LIMON)
library(here)
library(lme4)
library(Matrix)
library(tscount)
library(patchwork)
library(MASS)
library(matrixcalc)
library(gridExtra)
library(devtools)
library(miaSim)
library(reshape2)
library(corpcor)
2. Simulate Counts
Simulate a different starting covariance matrix and counts per
subject. ie each one has their own seed
# Initalize an empty list
####################################################
count_tables1 <- list()
# The Loop
####################################################
# Set the indices
for (i in 1:50) {
# Set the seed for each subject
set.seed(12345 + i)
# Create Starting covariance matrix
n <- 100^2
lower_values <- runif(n / 2, min = -0.7, max = -0.2)
upper_values <- runif(n / 2, min = 0.2, max = 0.7)
user_interactions <- c(lower_values, upper_values)
user_interactions <- sample(user_interactions)
user_A <- randomA(
n_species = 100,
interactions = user_interactions,
diagonal = -1.0,
connectance = 0.5,
scale_off_diagonal = 0.2,
symmetric = TRUE)
# Heatmap of the covariance matrix
color_palette <- colorRampPalette(c("red", "white", "blue"))(20)
heatmap(
user_A,
Colv = NA,
Rowv = NA,
col = color_palette,
zlim = c(-1, 1),
main = paste("Subject", i)
)
# Run the simulation
tse_glv <- simulateGLV(n_species = 100,
A = user_A,
t_start = 0,
t_store = 4,
stochastic = FALSE,
norm = FALSE,
error_variance = 0.1)
# Get the count table
sim_data <- tse_glv@assays@data@listData[["counts"]]
# Store the count table in the list
count_tables1[[i]] <- t(sim_data)
}


















































# Merge together
####################################################
# Combine the count tables
combined_count_table <- do.call(rbind, count_tables1)
# Add Rownames
rownames(combined_count_table) <- paste0("Sbj", rep(1:50,
each = nrow(count_tables1[[1]])), "_Time", 1:4)
3. Create fake metadata
Make Metadata and merge with the count data
# Set seed
set.seed(12345)
# Df 1 is Metadata
########################################################
meta_data <- expand.grid(Time = 1:4,ID = 1:50)
rownames(meta_data) <- rownames(combined_count_table)
# Df 2 is Metadata merged with Counts
########################################################
#Round off and increase
combined_count_table <- as.data.frame(combined_count_table + abs(min(combined_count_table)))
combined_count_table <- (combined_count_table)*10
meta_counts <- base::merge(meta_data, combined_count_table, by ="row.names", all = TRUE)
meta_counts <- column_to_rownames(meta_counts, "Row.names")
Plot the data to look at the differences among subjects
Graphs to Check
# Individual Species Plots
########################################################
# Pivot to long data
count_long <- tidyr::pivot_longer(meta_counts, cols = starts_with("sp"), names_to = "Species")
# Plot the data
count_long %>%
ggplot(aes(x = Time, y = value, colour = as.factor(ID),
group = as.factor(ID), linetype = as.factor(ID))) +
geom_line() +
geom_point() +
geom_jitter() +
ylab("Count") +
labs(linetype = "ID", color = "ID") +
facet_wrap(~ Species) + # Create a panel for each species
theme(legend.position = "none") +
ggtitle("Time Series Data")

# Plot only two
########################################################
# Pivot to long data
metacounts_filt <- meta_counts[,c(1:6)]
count_long <- tidyr::pivot_longer(metacounts_filt, cols = starts_with("sp"), names_to = "Species")
# Plot the data
count_long %>%
ggplot(aes(x = Time, y = value, colour = as.factor(ID),
group = as.factor(ID), linetype = as.factor(ID))) +
geom_line() +
geom_point() +
geom_jitter() +
ylab("Count") +
labs(linetype = "ID", color = "ID") +
facet_wrap(~ Species) + # Create a panel for each species
theme_linedraw() +
theme(axis.text.x = element_text(family = "arial",color = "black", size = 14),
axis.text.y = element_text(family = "arial",color = "black", size = 14),
axis.title.x = element_text(family = "arial",color = "black", size = 14),
axis.title.y = element_text(family = "arial",color = "black", size = 14),
legend.position = "none",
strip.text = element_text(face = "bold", family = "arial", color = "white", size = 15))

# Distribution of counts
########################################################
hist(as.matrix(combined_count_table),
breaks = 100, main = "Distribution of GLV Data", xlab = "Counts")

4. Save the data
Save the Counts
# Make more like counts
meta_counts[,3:102] <- round(meta_counts[,3:102])
write.csv(meta_counts, here("Data","GLV_SimData", "Dataset_2f", "GLV_100sp.csv"))
# Save just timepoint 2
meta_counts_filt <- meta_counts %>% filter(Time == 2)
write.csv(meta_counts_filt, here("Data","GLV_SimData", "Dataset_2f", "GLV_100sp_T2.csv"))
LS0tCnRpdGxlOiAiRGF0YXNldCAyZiAxMDAgU2ltdWxhdGlvbiIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKU2FtZSBhcyBEYXRhc2V0IDJmIGJ1dCBzYW1wbGUgc2l6ZSBvZiA1MCBhbmQgbnVtYmVyIG9mIHRheGEgdG8gMTAwLiBBbGwgZGlmZmVyZW50IGludGVyYWN0aW9uIG1hdHJpY2VzCgojIDEuIExvYWQgTGlicmFyeQoqKioKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGlncmFwaCkKbGlicmFyeShOQlpJTU0pCmxpYnJhcnkoU3BpZWNFYXNpKQpsaWJyYXJ5KExJTU9OKQpsaWJyYXJ5KGhlcmUpCmxpYnJhcnkobG1lNCkKbGlicmFyeShNYXRyaXgpCmxpYnJhcnkodHNjb3VudCkKbGlicmFyeShwYXRjaHdvcmspCmxpYnJhcnkoTUFTUykKbGlicmFyeShtYXRyaXhjYWxjKQpsaWJyYXJ5KGdyaWRFeHRyYSkKbGlicmFyeShkZXZ0b29scykKbGlicmFyeShtaWFTaW0pCmxpYnJhcnkocmVzaGFwZTIpCmxpYnJhcnkoY29ycGNvcikKYGBgCgoKCiMgMi4gU2ltdWxhdGUgQ291bnRzCgpTaW11bGF0ZSBhIGRpZmZlcmVudCBzdGFydGluZyBjb3ZhcmlhbmNlIG1hdHJpeCBhbmQgY291bnRzIHBlciBzdWJqZWN0LiBpZSBlYWNoIG9uZSBoYXMgdGhlaXIgb3duIHNlZWQKYGBge3J9CiMgSW5pdGFsaXplIGFuIGVtcHR5IGxpc3QKIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwpjb3VudF90YWJsZXMxIDwtIGxpc3QoKQoKIyBUaGUgTG9vcAojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCgojIFNldCB0aGUgaW5kaWNlcwpmb3IgKGkgaW4gMTo1MCkgewogIAogICMgU2V0IHRoZSBzZWVkIGZvciBlYWNoIHN1YmplY3QKICBzZXQuc2VlZCgxMjM0NSArIGkpICAKICAKICAjIENyZWF0ZSBTdGFydGluZyBjb3ZhcmlhbmNlIG1hdHJpeAogIG4gPC0gMTAwXjIgIAogIGxvd2VyX3ZhbHVlcyA8LSBydW5pZihuIC8gMiwgbWluID0gLTAuNywgbWF4ID0gLTAuMikgIAogIHVwcGVyX3ZhbHVlcyA8LSBydW5pZihuIC8gMiwgbWluID0gMC4yLCBtYXggPSAwLjcpCiAgdXNlcl9pbnRlcmFjdGlvbnMgPC0gYyhsb3dlcl92YWx1ZXMsIHVwcGVyX3ZhbHVlcykKICB1c2VyX2ludGVyYWN0aW9ucyA8LSBzYW1wbGUodXNlcl9pbnRlcmFjdGlvbnMpCiAgdXNlcl9BICA8LSByYW5kb21BKAogICAgbl9zcGVjaWVzID0gMTAwLAogICAgaW50ZXJhY3Rpb25zID0gdXNlcl9pbnRlcmFjdGlvbnMsCiAgICBkaWFnb25hbCA9IC0xLjAsCiAgICBjb25uZWN0YW5jZSA9IDAuNSwKICAgIHNjYWxlX29mZl9kaWFnb25hbCA9IDAuMiwKICAgIHN5bW1ldHJpYyA9IFRSVUUpCiAgCiAgIyBIZWF0bWFwIG9mIHRoZSBjb3ZhcmlhbmNlIG1hdHJpeAogIGNvbG9yX3BhbGV0dGUgPC0gY29sb3JSYW1wUGFsZXR0ZShjKCJyZWQiLCAid2hpdGUiLCAiYmx1ZSIpKSgyMCkKICBoZWF0bWFwKAogICAgdXNlcl9BLAogICAgQ29sdiA9IE5BLAogICAgUm93diA9IE5BLAogICAgY29sID0gY29sb3JfcGFsZXR0ZSwKICAgIHpsaW0gPSBjKC0xLCAxKSwgCiAgICBtYWluID0gcGFzdGUoIlN1YmplY3QiLCBpKQogICkKICAKICAjIFJ1biB0aGUgc2ltdWxhdGlvbgogIHRzZV9nbHYgPC0gc2ltdWxhdGVHTFYobl9zcGVjaWVzID0gMTAwLAogICAgICAgICAgICAgICAgICAgICAgICAgQSA9IHVzZXJfQSwKICAgICAgICAgICAgICAgICAgICAgICAgIHRfc3RhcnQgPSAwLCAKICAgICAgICAgICAgICAgICAgICAgICAgIHRfc3RvcmUgPSA0LAogICAgICAgICAgICAgICAgICAgICAgICAgc3RvY2hhc3RpYyA9IEZBTFNFLAogICAgICAgICAgICAgICAgICAgICAgICAgbm9ybSA9IEZBTFNFLAogICAgICAgICAgICAgICAgICAgICAgICAgZXJyb3JfdmFyaWFuY2UgPSAwLjEpCgogICMgR2V0IHRoZSBjb3VudCB0YWJsZQogIHNpbV9kYXRhIDwtIHRzZV9nbHZAYXNzYXlzQGRhdGFAbGlzdERhdGFbWyJjb3VudHMiXV0KICAKICAjIFN0b3JlIHRoZSBjb3VudCB0YWJsZSBpbiB0aGUgbGlzdAogIGNvdW50X3RhYmxlczFbW2ldXSA8LSB0KHNpbV9kYXRhKQp9CgojIE1lcmdlIHRvZ2V0aGVyCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMKCiMgQ29tYmluZSB0aGUgY291bnQgdGFibGVzCmNvbWJpbmVkX2NvdW50X3RhYmxlIDwtIGRvLmNhbGwocmJpbmQsIGNvdW50X3RhYmxlczEpCgojIEFkZCBSb3duYW1lcwpyb3duYW1lcyhjb21iaW5lZF9jb3VudF90YWJsZSkgPC0gcGFzdGUwKCJTYmoiLCByZXAoMTo1MCwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBlYWNoID0gbnJvdyhjb3VudF90YWJsZXMxW1sxXV0pKSwgIl9UaW1lIiwgMTo0KQpgYGAKCgoKCiMgMy4gQ3JlYXRlIGZha2UgbWV0YWRhdGEgIAoKCk1ha2UgTWV0YWRhdGEgYW5kIG1lcmdlIHdpdGggdGhlIGNvdW50IGRhdGEKYGBge3J9CiMgU2V0IHNlZWQKc2V0LnNlZWQoMTIzNDUpCiMgRGYgMSBpcyBNZXRhZGF0YQojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwptZXRhX2RhdGEgPC0gIGV4cGFuZC5ncmlkKFRpbWUgPSAxOjQsSUQgPSAxOjUwKQpyb3duYW1lcyhtZXRhX2RhdGEpIDwtIHJvd25hbWVzKGNvbWJpbmVkX2NvdW50X3RhYmxlKQoKCgojIERmIDIgaXMgTWV0YWRhdGEgbWVyZ2VkIHdpdGggQ291bnRzCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCiNSb3VuZCBvZmYgYW5kIGluY3JlYXNlCmNvbWJpbmVkX2NvdW50X3RhYmxlIDwtIGFzLmRhdGEuZnJhbWUoY29tYmluZWRfY291bnRfdGFibGUgKyBhYnMobWluKGNvbWJpbmVkX2NvdW50X3RhYmxlKSkpCmNvbWJpbmVkX2NvdW50X3RhYmxlIDwtIChjb21iaW5lZF9jb3VudF90YWJsZSkqMTAKbWV0YV9jb3VudHMgPC0gYmFzZTo6bWVyZ2UobWV0YV9kYXRhLCBjb21iaW5lZF9jb3VudF90YWJsZSwgYnkgPSJyb3cubmFtZXMiLCBhbGwgPSBUUlVFKQptZXRhX2NvdW50cyA8LSBjb2x1bW5fdG9fcm93bmFtZXMobWV0YV9jb3VudHMsICJSb3cubmFtZXMiKQpgYGAKCgpQbG90IHRoZSBkYXRhIHRvIGxvb2sgYXQgdGhlIGRpZmZlcmVuY2VzIGFtb25nIHN1YmplY3RzCgpHcmFwaHMgdG8gQ2hlY2sKYGBge3J9CiMgSW5kaXZpZHVhbCBTcGVjaWVzIFBsb3RzCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCiMgUGl2b3QgdG8gbG9uZyBkYXRhCmNvdW50X2xvbmcgPC0gdGlkeXI6OnBpdm90X2xvbmdlcihtZXRhX2NvdW50cywgY29scyA9IHN0YXJ0c193aXRoKCJzcCIpLCBuYW1lc190byA9ICJTcGVjaWVzIikKCiMgUGxvdCB0aGUgZGF0YQpjb3VudF9sb25nICU+JQogIGdncGxvdChhZXMoeCA9IFRpbWUsIHkgPSB2YWx1ZSwgY29sb3VyID0gYXMuZmFjdG9yKElEKSwKICAgICAgICAgICAgIGdyb3VwID0gYXMuZmFjdG9yKElEKSwgbGluZXR5cGUgPSBhcy5mYWN0b3IoSUQpKSkgKwogIGdlb21fbGluZSgpICsgCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2ppdHRlcigpICsKICB5bGFiKCJDb3VudCIpICsKICBsYWJzKGxpbmV0eXBlID0gIklEIiwgY29sb3IgPSAiSUQiKSArCiAgZmFjZXRfd3JhcCh+IFNwZWNpZXMpICsgICMgQ3JlYXRlIGEgcGFuZWwgZm9yIGVhY2ggc3BlY2llcwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKwogIGdndGl0bGUoIlRpbWUgU2VyaWVzIERhdGEiKQoKCiMgUGxvdCBvbmx5IHR3bwojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwojIFBpdm90IHRvIGxvbmcgZGF0YQptZXRhY291bnRzX2ZpbHQgPC0gbWV0YV9jb3VudHNbLGMoMTo2KV0KY291bnRfbG9uZyA8LSB0aWR5cjo6cGl2b3RfbG9uZ2VyKG1ldGFjb3VudHNfZmlsdCwgY29scyA9IHN0YXJ0c193aXRoKCJzcCIpLCBuYW1lc190byA9ICJTcGVjaWVzIikKCiMgUGxvdCB0aGUgZGF0YQpjb3VudF9sb25nICU+JQogIGdncGxvdChhZXMoeCA9IFRpbWUsIHkgPSB2YWx1ZSwgY29sb3VyID0gYXMuZmFjdG9yKElEKSwKICAgICAgICAgICAgIGdyb3VwID0gYXMuZmFjdG9yKElEKSwgbGluZXR5cGUgPSBhcy5mYWN0b3IoSUQpKSkgKwogIGdlb21fbGluZSgpICsgCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2ppdHRlcigpICsKICB5bGFiKCJDb3VudCIpICsKICBsYWJzKGxpbmV0eXBlID0gIklEIiwgY29sb3IgPSAiSUQiKSArCiAgZmFjZXRfd3JhcCh+IFNwZWNpZXMpICsgICMgQ3JlYXRlIGEgcGFuZWwgZm9yIGVhY2ggc3BlY2llcwogIHRoZW1lX2xpbmVkcmF3KCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGZhbWlseSA9ICJhcmlhbCIsY29sb3IgPSAiYmxhY2siLCBzaXplID0gMTQpLAogICAgICAgIGF4aXMudGV4dC55ID0gZWxlbWVudF90ZXh0KGZhbWlseSA9ICJhcmlhbCIsY29sb3IgPSAiYmxhY2siLCBzaXplID0gMTQpLAogICAgICAgIGF4aXMudGl0bGUueCA9IGVsZW1lbnRfdGV4dChmYW1pbHkgPSAiYXJpYWwiLGNvbG9yID0gImJsYWNrIiwgc2l6ZSA9IDE0KSwKICAgICAgICBheGlzLnRpdGxlLnkgPSBlbGVtZW50X3RleHQoZmFtaWx5ID0gImFyaWFsIixjb2xvciA9ICJibGFjayIsIHNpemUgPSAxNCksCiAgICAgICAgbGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiLAogICAgICAgIHN0cmlwLnRleHQgPSBlbGVtZW50X3RleHQoZmFjZSA9ICJib2xkIiwgZmFtaWx5ID0gImFyaWFsIiwgY29sb3IgPSAid2hpdGUiLCBzaXplID0gMTUpKQoKCiMgRGlzdHJpYnV0aW9uIG9mIGNvdW50cwojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwpoaXN0KGFzLm1hdHJpeChjb21iaW5lZF9jb3VudF90YWJsZSksIAogICAgIGJyZWFrcyA9IDEwMCwgbWFpbiA9ICJEaXN0cmlidXRpb24gb2YgR0xWIERhdGEiLCB4bGFiID0gIkNvdW50cyIpCmBgYAoKCiMgNC4gU2F2ZSB0aGUgZGF0YQoKU2F2ZSB0aGUgQ291bnRzCmBgYHtyfQojIE1ha2UgbW9yZSBsaWtlIGNvdW50cwptZXRhX2NvdW50c1ssMzoxMDJdIDwtIHJvdW5kKG1ldGFfY291bnRzWywzOjEwMl0pCndyaXRlLmNzdihtZXRhX2NvdW50cywgaGVyZSgiRGF0YSIsIkdMVl9TaW1EYXRhIiwgIkRhdGFzZXRfMmYiLCAiR0xWXzEwMHNwLmNzdiIpKQoKIyBTYXZlIGp1c3QgdGltZXBvaW50IDIKbWV0YV9jb3VudHNfZmlsdCA8LSBtZXRhX2NvdW50cyAlPiUgZmlsdGVyKFRpbWUgPT0gMikKd3JpdGUuY3N2KG1ldGFfY291bnRzX2ZpbHQsIGhlcmUoIkRhdGEiLCJHTFZfU2ltRGF0YSIsICJEYXRhc2V0XzJmIiwgIkdMVl8xMDBzcF9UMi5jc3YiKSkKYGBgCgoKCgo=